A Robust Learning Approach for Regression Models Based on Distributionally Robust Optimization

We present a Distributionally Robust Optimization (DRO) approach to estimate a robustified regression plane in a linear regression setting, when the observed samples are potentially contaminated with adversarially corrupted outliers. Our approach mitigates the impact of outliers by hedging against a family of probability distributions on the observed data, some of which assign very low probabilities to the outliers. The set of distributions under consideration are close to the empirical distribution in the sense of the Wasserstein metric. We show that this DRO formulation can be relaxed to a convex optimization problem which encompasses a class of models. By selecting proper norm spaces for the Wasserstein metric, we are able to recover several commonly used regularized regression models. We provide new insights into the regularization term and give guidance on the selection of the regularization coefficient from the standpoint of a confidence region. We establish two types of performance guarantees for the solution to our formulation under mild conditions. One is related to its out-of-sample behavior (prediction bias), and the other concerns the discrepancy between the estimated and true regression planes (estimation bias). Extensive numerical results demonstrate the superiority of our approach to a host of regression models, in terms of the prediction and estimation accuracies. We also consider the application of our robust learning procedure to outlier detection, and show that our approach achieves a much higher AUC (Area Under the ROC Curve) than M-estimation (Huber, 1964, 1973).

Given samples (x i , y i ),i = 1, … , N, we are interested in estimating β*. The Ordinary Least Squares (OLS) minimizes the sum of squared residuals ∑ i = 1 N y i − x i ′β 2 , and works well if all the N samples are generated from the underlying true model. However, when faced with adversarial perturbations in the training data, the OLS estimator will deviate from the true regression plane to reduce large residuals. Alternatively, one can choose to minimize the sum of absolute residuals ∑ i = 1 N |y i − x i ′β|, as done in Least Absolute Deviation (LAD), to mitigate the influence of large residuals. Another commonly used approach for hedging against outliers is M-estimation (Huber, 1964(Huber, , 1973, which minimizes a symmetric loss function ρ(·) of the residuals in the form ∑ i = 1 N ρ y i − x i ′β , which downweights the influence of samples with large absolute residuals. Several choices for ρ(·) include the Huber function (Huber, 1964(Huber, , 1973, the Tukey's Biweight function (Rousseeuw and Leroy, 2005), the logistic function (Coleman et al., 1980), the Talwar function (Hinich and Talwar, 1975), and the Fair function (Fair, 1974).
Both LAD and M-estimation are not resistant to large deviations in the predictors. For contamination present in the predictor space, high breakdown value methods are required. Examples include the Least Median of Squares (LMS) , which minimizes the median of the absolute residuals, the Least Trimmed Squares (LTS) (Rousseeuw, 1985), which minimizes the sum of the q smallest squared residuals, and S-estimation (Rousseeuw and Yohai, 1984), which has a higher statistical efficiency than LTS with the same breakdown value. A combination of the high breakdown value method and M-estimation is the MM-estimation (Yohai, 1987). It has a higher statistical efficiency than S-estimation. We refer the reader to the book of Rousseeuw and Leroy (2005) for a detailed description of these robust regression methods.
The aforementioned robust estimation procedures focus on modifying the objective function in a heuristic way with the intent of minimizing the effect of outliers. A more rigorous line of research explores the underlying stochastic program that leads to the sample-based estimation procedures. For example, the OLS objective can be viewed as minimizing the expected squared residual under the uniform empirical distribution over the samples. It has been well recognized that optimizing under the empirical distribution yields estimators that are sensitive to perturbations in the data and suffer from overfitting. The reason is that when the data (x, y) are adversarially corrupted by outliers, the observed samples do not represent well the true underlying distribution of the data. Yet, the samples are typically the only information available. Instead of equally weighting all the samples as in the empirical distribution, we may wish to include more informative distributions that "drive out" the corrupted samples. One way to realize this is to hedge the expected loss against a family notion of closeness between two points, which in the context of image retrieval represents the perceptual similarity in color.
Our DRO problem minimizes the worst-case absolute residual over a Wasserstein ball of distributions, and could be relaxed to the following form: where ϵ is the radius of the Wasserstein ball, and ∥ · ∥ * is the dual norm of the norm space where the Wasserstein metric is defined on. Formulation (1) incorporates a wide class of models whose specific form depends on the notion of transportation cost embedded in the Wasserstein metric (see Section 2). Although the Wasserstein DRO formulation simply reduces to regularized regression models, we want to emphasize a few new insights brought by this methodology. First, the regularization term controls the conservativeness of the Wasserstein set, or the amount of ambiguity in the data, which differentiates itself from the heuristically added regularizers in traditional regression models that serve the purpose of preventing overfitting, error/variance reduction, or sparsity recovery. Second, the regularization term is determined by the dual norm of the regression coefficient, which controls the growth rate of the ℓ 1 -loss function, and the radius of the Wasserstein set. This connection provides guidance on the selection of the regularization coefficient and may lead to significant computational savings compared to cross-validation. DRO essentially enables new and more accurate interpretations of the regularizer, and establishes its dependence on the growth rate of the loss, the underlying metric space and the reliability of the observed samples.
The connection between robustness and regularization has been established in several works. The earliest one may be credited to El Ghaoui and Lebret (1997), who show that minimizing the worst-case squared residual within a Frobenius norm-based perturbation set is equivalent to Tikhonov regularization. In more recent works, using properly selected uncertainty sets, Xu et al. (2010)  Our work is motivated by the problem of identifying patients who receive an abnormally high radiation exposure in CT exams, given the patient characteristics and exam-related variables (Chen et al., 2018). This could be casted as an outlier detection problem; specifically, estimating a robustified regression plane that is immunized against outliers and learns the underlying true relationship between radiation dose and the relevant predictors. We focus on robust learning of the parameter in regression models under distributional perturbations residing within a Wasserstein ball. While the applicability of the Wasserstein DRO methodology is not restricted to regression analysis (Sinha et al., 2017;Gao et al., 2017;Shafieezadeh-Abadeh et al., 2017), or a particular form of the loss function (as long as it satisfies certain smoothness conditions (Gao et al., 2017)), we focus on the absolute residual loss in linear regression in light of our motivating application and for the purpose of enhancing robustness. Our contributions can be summarized as follows:

1.
We develop a DRO approach to robustify linear regression using an ℓ 1 loss function and an ambiguity set around the empirical distribution of the training samples defined based on the Wasserstein metric. The formulation is general enough to include any norm-induced Wasserstein metric and incorporate additional regularization constraints on the regression coefficients (e.g., ℓ 1 -norm constraints). It provides an intuitive connection between the amount of ambiguity allowed and a regularization penalty term in the robust formulation, which provides a natural way to adjust the latter.

2.
We establish novel performance guarantees on both the out-of-sample loss (predictionbias) and the discrepancy between the estimated and the true regression coefficients (estimation bias). Our guarantees elucidate the role of the regularizer, which is related to the dual norm of the regression coefficients, in bounding the biases and are in concert with the theoretical foundation that leads to the regularized problem. The generalization error bound, in particular, builds a connection between the loss function and the form of the regularizer via Rademacher complexity, providing a rigorous explanation for the commonly observed good out-of-sample performance of regularized regression. On the other hand, the estimation error bound corroborates the validity of the ℓ 1 -loss function, which tends to incur a lower estimation bias than other candidates such as the ℓ 2 and ℓ ∞ losses. Our results are novel in the robust regression setting and different from earlier work in the DRO literature, enabling new perspectives and interpretations of the norm-based regularization, and providing justifications for the ℓ 1 -loss-based learning algorithms.

3.
We empirically explore three important aspects of the Wasserstein DRO formulation, including the advantages of the ℓ 1 -loss function, the selection of a proper norm for the Wasserstein metric, and the implication of penalizing the extended regression coefficient (−β, 1), by comparing with a series of regression models on a number of synthetic datasets. We show the superiority of the Wasserstein DRO approach, presenting a thorough analysis under four different experimental setups. We also consider the application of our methodology to outlier detection and compare with M-estimation in terms of the ability of identifying outliers (ROC (Receiver Operating Characteristic) curves). The Wasserstein DRO formulation achieves significantly higher AUC (Area Under Curve) values.
The rest of the paper is organized as follows. In Section 2, we introduce the Wasserstein metric and derive the general Wasserstein DRO formulation in a linear regression framework. Section 3 establishes performance guarantees for both the general formulation and the special case where the Wasserstein metric is defined on the ℓ 1 -norm space.
Numerical experimental results are presented in Section 4. We conclude the paper in Section 5.

Notational conventions:
We use boldfaced lowercase letters to denote vectors, ordinary lowercase letters to denote scalars, boldfaced uppercase letters to denote matrices, and calligraphic capital letters to denote sets. E denotes expectation and ℙ probability of an event. All vectors are column vectors. For space saving reasons, we write x = (x 1 , … , x dim(x) ) to denote the column vector We use prime to denote the transpose of a vector, ∥·∥ for the general norm operator, ∥·∥ 2 for the ℓ 2 norm, ∥ · ∥ 1 for the ℓ 1 norm, and ∥ · ∥ ∞ for the infinity norm. P(Z) denotes the set of probability measures supported on Z. e i denotes the i-th unit vector, e the vector of ones, 0 a vector of zeros, and I the identity matrix. Given a norm ∥ · ∥ on ℝ m , the dual norm ∥ · ∥ * is defined as: ∥θ∥ * ≜ sup ∥z∥≤1 θ′z. For a function h(z), its convex conjugate h*(·) is defined as: h*(θ) ≜ sup z∈dom h {θ′z − h(z)}, where dom h denotes the domain of the function h.

Problem Statement and Justification of Our Formulation
Consider a linear regression problem where we are given a predictor/feature vector x ∈ ℝ m − 1 , and a response variable y ∈ ℝ. Our goal is to obtain an accurate estimate of the regression plane that is robust with respect to the adversarial perturbations in the data. We consider an ℓ 1 -loss function h β (x, y) ≜ |y − x′β|, motivated by the observation that the absolute loss function is more robust to large residuals than the squared loss (see Fig. 1). Moreover, the estimation error analysis presented in Section 3.2 suggests that the ℓ 1 -loss function leads to a smaller estimation bias than others. Our Wasserstein DRO problem using the ℓ 1 -loss function is formulated as: where β is the regression coefficient vector that belongs to some set ℬ. ℬ could be ℝ m − 1 , or ℬ = β: β 1 ≤ l if we wish to induce sparsity, with l being some pre-specified number. ℚ is the probability distribution of (x, y), belonging to some set Ω which is defined as: W p ℚ, ℙ N is the order-p Wasserstein distance between ℚ and ℙ N (see definition in (3)), with ℙ N the uniform empirical distribution over samples. The formulation in (2) is robust since it minimizes over the regression coefficients the worst case expected loss, that is, the expected loss maximized over all probability distributions in the ambiguity set Ω.
Before deriving a tractable reformulation for (2), let us first define the Wasserstein metric. Let (Z, s) be a metric space where Z is a set and s is a metric on Z. The Wasserstein metric of order p ≥ 1 defines the distance between two probability distributions ℚ 1 and ℚ 2 in the following way: W p ℚ 1 , ℚ 2 ≜ min Π ∈ P(Z × Z) ∫ Z × Z s x 1 , y 1 , x 2 , y 2 p Π d x 1 , y 1 , d x 2 , y 2 where Π is the joint distribution of (x 1 , y 1 ) and (x 2 , y 2 ) with marginals ℚ 1 and ℚ 2 , respectively. The Wasserstein distance between ℚ 1 and ℚ 2 represents the cost of an optimal mass transportation plan, where the cost is measured through the metric s. The order p should be selected in such a way as to ensure that the worst-case expected loss is meaningfully defined, i.e., Notice that the ambiguity set Ω is centered at the empirical distribution ℙ N and has radius ϵ.
A formal definition of the growth rate is due to Gao and Kleywegt (2016), which takes the limit of (7) as s((x 1 , y 1 ),(x 2 , y 2 )) → ∞, to eliminate its dependence on (x, y). One important aspect they have pointed out is that when the growth rate of the loss function is infinite, strong duality for the worst-case problem sup ℚ ∈ Ω E ℚ ℎ β (x, y) fails to hold, in which case the DRO problem (2) becomes intractable. Assuming that the metric s is induced by some norm ∥ · ∥, the bounded growth rate requirement is expressed as follows: limsup ‖ x 1 , y 1 − x 2 , y 2 ‖ ∞ ℎ β x 1 , y 1 − ℎ β x 2 , y 2 x 1 , y 1 − x 2 , y 2 p ≤ limsup ‖ x 1 , y 1 − x 2 , y 2 ‖ ∞ y 1 − x 1 ′ β − y 2 − x 2 ′ β x 1 , y 1 − x 2 , y 2 p ≤ lim sup ‖ x 1 , y 1 − x 2 , y 2 ‖ ∞ x 1 , y 1 − x 2 , y 2 ( − β, 1) * x 1 , y 1 − x 2 , y 2 p < ∞, where ∥·∥ * is the dual norm of ∥·∥, and the second inequality is due to the Cauchy-Schwarz inequality. Notice that by taking p = 1, (8) is equivalently translated into the condition that ∥(−β, 1)∥ * < ∞, which we will see in Section 3 is an essential requirement to guarantee a good generalization performance for the Wasserstein DRO estimator. The growth rate essentially reveals the underlying metric space used by the Wasserstein distance. Taking p > 1 leads to zero growth rate in the limit of (8), which is not desirable since it removes the Wasserstein ball structure from our formulation and renders it an optimization problem over a singleton distribution. This will be made more clear in the following analysis. We thus choose the order-1 Wasserstein metric with s being induced by some norm ∥ · ∥ to define our DRO problem.
Next, we will discuss how to convert (2) into a tractable formulation. Suppose we have N independently and identically distributed realizations of (x, y), denoted by (x i , y i ), i = 1, … , N. We make the assumption that (x, y) comes from a mixture of two distributions, with probability q from the outlying distribution ℙ out and with probability 1 − q from the true distribution ℙ. Recall that ℙ N is the discrete uniform distribution over the N samples. Our goal is to generate estimators that are consistent with the true distribution ℙ. We claim that when q is small, if the Wasserstein ball radius ϵ is chosen judiciously, the true distribution ℙ will be included in the set Ω while the outlying distribution ℙ out will be excluded. To see this, consider a simple example where ℙ is a discrete distribution that assigns equal probability to 10 data points equally spaced between 0.1 and 1, and ℙ out assigns probability robust to the adversarial perturbations. Moreover, as q becomes smaller, the gap between the red and blue lines becomes larger. One implication from this observation is that as the data becomes purer, the radius of the Wasserstein ball tends to be smaller, and the confidence in the observed samples is higher. For large q values, the DRO formulation seems to fail.
However, as outliers are defined to be the data points that do not conform to the majority of data, we can safely claim that ℙ out is the distribution of the minority and q is always below 0.5.
Through (9), we can relax problem (2) by minimizing the right hand side of (9) instead of the worst-case expected loss. Moreover, as shown in Esfahani and Kuhn (2015), (9) becomes an equality when Z = ℝ m . In Theorem 2.1, we compute the value of κ(β) for the specific ℓ 1 loss function we use. The proof of this Theorem and all results hereafter are included in Appendix A.
Due to Theorem 2.1, (2) could be formulated as the following optimization problem: Note that the regularization term of (10) is the product of the growth rate of the loss and the Wasserstein ball radius. The growth rate is closely related to the way the Wasserstein metric defines the transportation costs on the data (x, y). As mentioned earlier, a zero growth rate diminishes the effect of the Wasserstein distributional uncertainty set, and the resulting formulation would simply be an empirical loss minimization problem. The parameter ϵ controls the conservativeness of the formulation, whose selection depends on the sample size, the dimensionality of the data, and the confidence that the Wasserstein ball contains the true distribution (see eq. (8) in Esfahani and Kuhn, 2015). Roughly speaking, when the sample size is large enough, and for a fixed confidence level, ϵ is inversely proportional to N 1/m . polyhedron with convex quadratic inequalities, (10) is a convex quadratic problem which can be solved to optimality very efficiently. Specifically, it could be converted to: When the Wasserstein metric is defined using ∥·∥ 1 and the set ℬ is a polyhedron, (10) is a linear programming problem: More generally, when the coordinates of (x, y) differ from each other substantially, a properly chosen, positive definite weight matrix M ∈ ℝ m × m could scale correspondingly different coordinates of (x, y) by using the M-weighted norm: It can be shown that (10) in this case becomes: We note that this Wasserstein DRO framework could be applied to a broad class of loss functions and the tractable reformulations have been derived in Shafieezadeh-Abadeh et al. (2017); Gao et al. (2017) for regression and classification models. We adopt the absolute residual loss in this paper to enhance the robustness of the formulation, which is the focus of our work and serves the purpose of estimating robust parameters that are immunized against perturbations/outliers. Notice that (10) coincides with the regularized LAD models (Pollard, 1991;Wang et al., 2006), except that we are regularizing a variant of the regression coefficient. We would like to highlight several novel viewpoints that are brought by the Wasserstein DRO framework and justify the value and novelty of (10). First, (10) is obtained as an outcome of a fundamental DRO formulation, which enables new interpretations of the regularizer from the standpoint of distributional robustness, and provides rigorous theoretical foundation on why the ℓ 2 -regularizer prevents overfitting to the training data. The regularizer could be seen as a control over the amount of ambiguity in the data and reveals the reliability of the contaminated samples. Second, the geometry of the Wasserstein ball is embedded in the regularization term, which penalizes the regression coefficient on the dual Wasserstein space, with the magnitude of penalty being the radius of the ball. This offers an intuitive interpretation and provides guidance on how to set the regularization coefficient. Moreover, different from the traditional regularized LAD models that directly penalize the regression coefficient β, we regularize the vector (−β, 1), where the 1 takes into account the transportation cost along the y direction. Penalizing only on β corresponds to an infinite transportation cost along y. Our model is more general in this sense, and establishes the connection between the metric space on data and the form of the regularizer.

Performance Guarantees
Having obtained a tractable reformulation for the Wasserstein DRO problem, we next establish guarantees on the predictive power and estimation quality for the solution to (10). Two types of results will be presented in this section, one of which bounds the prediction bias of the estimator on new, future data (given in Section 3.1). The other one that bounds the discrepancy between the estimated and true regression planes (estimation bias), is given in Section 3.2.

Out-of-Sample Performance
In this subsection we investigate generalization characteristics of the solution to (10), which involves measuring the error generated by our estimator on a new random sample (x, y).
We would like to obtain estimates that not only explain the observed samples well, but, more importantly, possess strong generalization abilities. The derivation is mainly based on Rademacher complexity (see Bartlett and Mendelson, 2002), which is a measurement of the complexity of a class of functions. We would like to emphasize the applicability of such a proof technique to general loss functions, as long as their empirical Rademacher complexity could be bounded. The bound we derive for the prediction bias depends on both the sample average loss (the training error) and the dual norm of the regression coefficient (the regularizer), which corroborates the validity and necessity of our regularized formulation. Moreover, the generalization result also builds a connection between the loss function and the form of the regularizer via Rademacher complexity, which enables new insights into the regularization term and explains the commonly observed good out-of-sample performance of regularized regression in a rigorous way. We first make several mild assumptions that are needed for the generalization result.
Assumption A-The norm of the uncertainty parameter (x, y) is bounded above almost surely, i.e., ∥(x, y)∥ ≤ R.
Assumption B-The dual norm of (−β, 1) is bounded above within the feasible region, namely, Under these two assumptions, the absolute loss could be bounded via the Cauchy-Schwarz inequality.
With the above result, the idea is to bound the generalization error using the empirical Rademacher complexity of the following class of loss functions: We need to show that the empirical Rademacher complexity of ℋ, denoted by ℛ N (ℋ), is upper bounded. The following result, similar to Lemma 3 in Bertsimas et al. (2015), provides a bound that is inversely proportional to the square root of the sample size.
Let β be an optimal solution to (10), obtained using the samples (x i , y i ), i = 1, … , N. Suppose we draw a new i.i.d. sample (x, y). In Theorem 3.3 we establish bounds on the error y − x′β .

Theorem 3.3-Under
Assumptions A and B, for any 0 < δ < 1, with probability at least 1 − δ with respect to the sampling, and for any ζ > 2BR There are two probability measures in the statement of Theorem 3.3. One is related to the new data (x, y), while the other is related to the samples (x 1 , y 1 ), … , (x N , y N ). The expectation in (14) (and the probability in (15)) is taken w.r.t. the new data (x, y). For a given set of samples, (14) (and (15)) holds with probability at least 1 − δ w.r.t. the measure of samples. Theorem 3.3 essentially says that given typical samples, the expected loss on new data using our Wasserstein DRO estimator could be bounded above by the average sample loss plus extra terms that depend on the supremum of ∥(−β, 1) ∥ * (our regularizer), and are proportional to 1/ N. This result validates the dual norm-based regularized regression from the perspective of generalization ability, and could be generalized to any bounded loss function. It also provides implications on the form of the regularizer. For example, if given an ℓ 2 -loss function, the dependency on B for the generalization error bound will be of the form B 2 , which suggests using ( − β, 1) * 2 as a regularizer, reducing to a variant of ridge regression (Hoerl and Kennard, 1970) for ∥ · ∥ 2 induced Wasserstein metric.
We also note that the upper bounds in (14) and (15) do not depend on the dimension of (x, y). This dimensionality-free characteristic implies direct applicability of our Wasserstein approach to high-dimensional settings and is particularly useful in many real applications where, potentially, hundreds of features may be present. Theorem 3.3 also provides guidance on the number of samples that are needed to achieve satisfactory out-of-sample performance.
Corollary 3.4-Suppose β is the optimal solution to (10). For a fixed confidence level δ and some threshold parameter τ ≥ 0, to guarantee that the percentage difference between the expected absolute loss on new data and the sample average loss is less than τ, that is, the sample size N must satisfy Corollary 3.5-Suppose β is the optimal solution to (10). For a fixed confidence level δ, some τ ∈ (0, 1) and γ ≥ 0, to guarantee that the sample size N must satisfy provided that τ · γ + τ − 1 > 0.
In Corollaries 3.4 and 3.5, the sample size is inversely proportional to both δ and τ, which is reasonable since the more confident we want to be, the more samples we need. Moreover, the smaller τ is, the stricter a requirement we impose on the performance, and thus more samples are needed.

Discrepancy between Estimated and True Regression Planes
In addition to the generalization performance, we are also interested in the accuracy of the estimator. In this section we seek to bound the difference between the estimated and true regression coefficients, under a certain distributional assumption on (x, y). Throughout the section we will use β to denote the estimated regression coefficients, obtained as an optimal solution to (18), and β* for the true (unknown) regression coefficients. The bound we will derive turns out to be related to the Gaussian width (see definition in the Appendix) of the unit ball in ∥ · ∥ ∞ , the sub-Gaussian norm of the uncertainty parameter (x, y), as well as the geometric structure of the true regression coefficients. We note that this proof technique may be applied to several other loss functions, e.g., ℓ 2 and ℓ ∞ losses, with slight modifications. However, we will see that the ℓ 1 -loss function incurs a relatively low estimation bias compared to others, further demonstrating the superiority of our absolute error minimization formulation.
To facilitate the analysis, we will use the following equivalent form of problem (10): where Z = [(x 1 , y 1 ), … , (x N , y N )] is the matrix with columns (x i , y i ),i = 1, … , N, and γ N is some exogenous parameter related to ϵ. One can show that for properly chosen γ N , (18) produces the same solution with (10) (Bertsekas, 1999). (18) is similar to (11) in Chen and Banerjee (2016), with the difference lying in that we impose a constraint on the error instead of the gradient, and we consider a more general notion of norm on the coefficient. On the other hand, due to their similarity, we will follow the line of development in Chen and Banerjee (2016). Still, our analysis is self-contained and the bound we obtain is in a different form, which provides meaningful insights into our specific problem. We list below the assumptions that are needed to bound the estimation error.

Assumption D (Restricted Eigenvalue Condition)-For some set
A β* = cone v | −β*, 1 + v * ≤ −β*, 1 * ∩ S m and some positive scalar α, where S m is the unit sphere in the m-dimensional Euclidean space, where S m denotes the unit sphere in the m-dimensional Euclidean space.
Assumption F-(x, y) is a centered sub-Gaussian random vector (see definition in the Appendix), i.e., it has zero mean and satisfies the following condition: (x, y) ψ 2 = sup u ∈ S m (x, y)′u ψ 2 ≤ μ .
Notice that both α in Assumption D and γ N in Assumption E are related to the random observation matrix Z. A probabilistic description for these two quantities will be provided later. We next present a preliminary result, similar to Lemma 2 in Chen and Banerjee (2016), that bounds the ℓ 2 -norm of the estimation bias in terms of a quantity that is related to the geometric structure of the true coefficients. This result gives a rough idea on the factors that affect the estimation error, and shows the advantages of using the ℓ 1 -loss from the perspective of its dual norm. The bound derived in Theorem 3.6 is crude in the sense that it is a function of several random parameters that are related to the random observation matrix Z. This randomness will be described in a probabilistic way in the subsequent analysis.
Notice that the bound in (19) does not explicitly depend on the sample size N. If we change to the ℓ 2 -loss function, problem (18) will become: The proof of Theorem 3.6 still applies with slight modification. We will find out that in the case of ℓ 2 -loss, the estimation error bound takes the following form: Similarly, the ℓ ∞ -loss, which considers only the maximum absolute loss among the samples, turns (18) into: The corresponding bound becomes: We see that by using either ℓ 2 or ℓ ∞ -loss, an explicit dependency on N is introduced. As a result, the estimation error bounds become worse. The reason is that for the ℓ 1 -loss function, its dual norm operator only picks out the maximum absolute coordinate and thus avoids the dependence on the dimension, which in our case is the sample size (see Eq. (28)), whereas other norms, e.g., ℓ 2 -norm, sum over all the coordinates and thus introduce a dependence on N.
As mentioned earlier, (19) provides a random upper bound, revealed in α and γ N , that depends on the randomness in Z. We therefore would like to replace these two parameters by non-random quantities. The α acts as the minimum eigenvalue of the matrix ZZ′ restricted to a subspace of ℝ m , and thus a proper substitute should be related to the minimum eigenvalue of the covariance matrix of (x, y), i.e., the Γ matrix (cf. Assumption G), given that (x, y) is zero mean. See Lemmas 3.7, 3.8 and 3.9 for the derivation.
Lemma 3.7-Consider the set A Γ = w ∈ S m | Γ −1/2 w ∈ cone A β* , where A β* is defined as in Theorem 3.6, and Γ = E (x, y)(x, y)′ . Under Assumptions F and G, when the sample size N ≥ C 1 μ 4 w A Γ 2 , where μ = μ 1 λ min , and w A Γ is the Gaussian width of A Γ , with probability at least 1 − exp −C 2 N /μ 4 , we have v′ZZ′v ≥ N 2 v′Γv, ∀v ∈ A β* , where C 1 and C 2 are positive constants.
Note that the sample size requirement stated in Lemma 3.7 depends on the Gaussian width of A Γ , where A Γ relates to A β* . The following lemma shows that their Gaussian widths are also related. This relation is built upon the square root of the eigenvalues of Γ, which measures the extent to which A Γ expands A β* . (2016))-Let μ 0 be the ψ 2 -norm of a standard Gaussian random vector g ∈ ℝ m , and A Γ , A β* be defined as in Lemma 3.7. Then, under Assumption G, w A Γ ≤ C 3 μ 0 λ max λ min w A β* + 3 , for some positive constant C 3 .

Lemma 3.8 (Lemma 4 in Chen and Banerjee
Combining Lemmas 3.7 and 3.8, and expressing the covariance matrix Γ using its eigenvalues, we arrive at the following result.

Corollary 3.9-Under
Assumptions F and G, and the conditions in Lemmas 3.7 and 3.8, where C 1 and C 2 are positive constants.
Next, we derive the smallest possible value of γ N such that β* is feasible. The derivation uses the dual norm operator of the ℓ 1 -loss, resulting in a bound that depends on the Gaussian width of the unit ball in the dual norm space (∥ · ∥ ∞ ). See Lemma 3.10 for details. where ℬ u is the unit ball of norm ∥·∥ ∞ , ρ = sup v ∈ ℬ u v 2 , and C 4 , C 5 , C positive constants.
We note that for other loss functions, e.g., the ℓ 2 and ℓ ∞ losses, similar results can be obtained, where ℬ u is defined to be the unit ⋅ * loss -ball in ℝ m , with ⋅ * loss being the dual norm of the loss. Combining Theorem 3.6, Corollary 3.9 and Lemma 3.10, we have the following main performance guarantee result that bounds the estimation bias of the solution to (18).
From (20) we see that the bias is decreased as the sample size increases and the uncertainty embedded in (x, y) (revealed in R and μ) is reduced. The estimation error bound depends on the geometric structure of the true coefficients, defined using the dual norm space of the Wasserstein metric, the Gaussian width of the unit ⋅ * loss -ball in ℝ m , and the minimum eigenvalue of the covariance matrix of (x, y), with a convergence rate 1/N for the ℓ 1 -loss we applied. As mentioned earlier, other loss functions may incur a dependence on N in the numerator of the bound, thus resulting in a slower convergence rate, which substantiates the benefit of using an ℓ 1 -loss function.

Simulation Experiments on Synthetic Datasets
In this section we will explore the robustness of the Wasserstein formulation in terms of its Absolute Deviation (AD) loss function and the dual norm regularizer on the extended regression coefficient (−β, 1). Recall that our Wasserstein formulation is in the following form: We will focus on the following three aspects of this formulation:

1.
How to choose a proper norm ∥ · ∥ for the Wasserstein metric?

3.
What is the advantage of the AD loss compared to the Squared Residuals (SR) loss?
To answer Question 1, we will connect the choice of ∥ · ∥ for the Wasserstein metric with the characteristics/structures of the data (x, y). Specifically, we will design two sets of experiments, one with a dense regression coefficient β*, where all coordinates of x play a role in determining the value of the response y, and another with a sparse β* implying that only a few predictors are relevant/important in predicting y. Two Wasserstein formulations will be tested and compared, one induced by the ∥ · ∥ 2 (Wasserstein ℓ 2 ), which leads to an ℓ 2 -regularizer in (21), and the other one induced by the ∥·∥ ∞ (Wasserstein ℓ ∞ ) and resulting in an ℓ 1 -regularizer in (21). Intuitively, and based on the past experience in implementing the regularization techniques, the Wasserstein ℓ 2 should outperform the Wasserstein ℓ ∞ in the dense setting, while in the sparse setting, the reverse is true. Researchers have well identified the sparsity inducing property of the ℓ 1 -regularizer and provided a nice geometrical interpretation for it (Friedman et al., 2001). Here, we try to offer a different explanation from the perspective of the Wasserstein DRO formulation, through projecting the sparsity of β* onto the (x, y) space and establishing a sparse distance metric that only extracts a subset of coordinates from (x, y) to measure the closeness between samples.
For the second question, we first note that if the Wasserstein metric is induced by the following metric s c : which is the ℓ 1 -regularized LAD (see proof in the Appendix). It follows that regularizing over β implies an infinite transportation cost along y. In other words, for two data points (x 1 , y 1 ) and (x 2 , y 2 ), if y 1 ≠ y 2 , then they are considered to be infinitely far away. By contrast, our Wasserstein formulation, which regularizes over the extended regression coefficient (−β, 1), stems from a finite cost along y that is equally weighted with x. We will see the disadvantages of penalizing only β in the analysis of the experimental results.
To answer Question 3, we will compare against several commonly used regression models that employ the SR loss function, e.g., ridge regression (Hoerl and Kennard, 1970), LASSO (Tibshirani, 1996), and Elastic Net (EN) (Zou and Hastie, 2005). We will also compare against M-estimation (Huber, 1964(Huber, , 1973, which uses a variant of the SR loss and is equivalent to solving a weighted least squares problem, where the weights are determined by the residuals. These models will be compared under two different experimental setups, one involving adversarial perturbations in both x and y, and the other with perturbations only in x. The purpose is to investigate the behavior of these approaches when the noise in y is substantially reduced. As shown by Fig. 1, compared to the SR loss, the AD loss is less vulnerable to large residuals, and hence, it is advantageous in the scenarios where large perturbations appear in y. We are interested in studying whether its performance is consistently good when the corruptions appear mainly in x. We next describe the data generation process. Each training sample has a probability q of being drawn from the outlying distribution, and a probability 1 − q of being drawn from the true (clean) distribution. Given the true regression coefficient β*, we generate the training data as follows: • Generate a uniform random variable on [0, 1]. If it is no larger than 1 − q, generate a clean sample as follows:

1.
Draw the predictor x ∈ ℝ m − 1 from the normal distribution N m−1 (0, Σ), where Σ is the covariance matrix of x, which is just the top left block of the matrix Γ in Assumption G. Specifically, Γ = E (x, y)(x, y)′ is equal to Γ = Σ Σβ* β* ′Σ β* ′Σβ* + σ 2 , with σ 2 being the variance of the noise term. In our implementation, Σ has diagonal elements equal to 1 (unit variance) and off-diagonal elements equal to ρ, with ρ the correlation between predictors.
• Otherwise, depending on the experimental setup, generate an outlier that is either: -Abnormal in both x and y, with outlying distribution: 1.

•
Repeat the above procedure for N times, where N is the size of the training set.
To test the generalization ability of various formulations, we generate a test dataset containing M samples from the clean distribution. It is worth noting that only clean samples are included in the test set, since we only care about the prediction accuracy on clean data points, and our estimator is supposed to be consistent with the clean distribution and stay away from the outlying one. We are interested in studying the performance of various methods as the following factors are varied: which is equally spaced between 0.05 and 2 on a log scale.
The performance metrics we use include: • Mean Squared Error (MSE) on the test dataset, which is defined to be • Relative Risk (RR) of β defined as: • Relative Test Error (RTE) of β defined as: • Proportion of Variance Explained (PVE) of β defined as: For the metrics that evaluate the accuracy of the estimator, i.e., the RR, RTE and PVE, we list below two types of scores, one achieved by the best possible estimator β = β*, called the perfect score, and the other one achieved by the null estimator β = 0, called the null score.
• RR: a perfect score is 0 and the null score is 1.
• RTE: a perfect score is 1 and the null score is SNR+1.
• PVE: a perfect score is sNR sNR+1 , and the null score is 0.
During the training process, all the regularization parameters are tuned on a separate validation dataset. Specifically, we divide all the N training samples into two sets, dataset 1 and dataset 2 (validation set). For a pre-specified range of values for the penalty parameters, dataset 1 is used to train the models and derive β , and the performance of β is evaluated on dataset 2. We choose the regularization parameter that yields the minimum Median Absolute Deviation (MAD) on the validation set. Using MAD as a selection criterion serves to hedge against the potentially large noise in the validation samples. As to the range of values for the tuned parameters, we borrow ideas from Hastie et al. (2017), where the LASSO was tuned over 50 values ranging from λ m = ∥X′y∥ ∞ to a small fraction of λ m on a log scale, with X ∈ ℝ N × m − 1 the design matrix whose i-th row is X i ′, and y = (y 1 , … , y N ) the response vector. In our experiments, this range is properly adjusted for procedures that use the AD loss. Specifically, for Wasserstein ℓ 2 and ℓ ∞ , ℓ 1 -and ℓ 2 -regularized LAD, the range of values for the regularization parameter is: exp lin log 0.005 * X′y ∞ , log X′y ∞ , 50 , where lin(a, b, n) is a function that takes in scalars a, b and n (integer) and outputs a set of n values equally spaced between a and b; the exp function is applied elementwise to a vector.
The square root operator is in consideration of the AD loss that is the square root of the SR loss if evaluated on a single sample.
The regularization coefficient ϵ in formulation (10), which is the radius of the Wasserstein ball, allows for a more efficient tuning procedure. It has been noted in Esfahani and Kuhn (2015) that for a large enough sample size, ϵ is inversely proportional to N 1/m . This proportionality could be used as a guidance on setting ϵ, where only the proportional factor needs to be tuned (using cross-validation or a separate validation dataset as described earlier). In our implementation, given the small size of the simulated datasets, we will still adopt the validation dataset approach to tune the regularization parameter.

Dense β*, outliers in both x and y
In this subsection, we choose a dense regression coefficient β*, set the intercept β 0 * = 0.3, and the coefficient for each predictor x i to be β i * = 0.5, i = 1, … , 20. The adversarial perturbations are present in both x and y. Specifically, the outlying distribution is described by: 1.
We generate 10 datasets consisting of N = 100, M = 60 observations. The probability of a training sample being drawn from the outlying distribution is q = 30%. The mean values of the performance metrics (averaged over the 10 datasets), as we vary the SNR and the correlation between predictors, are shown in Figs. 3 and 4. Note that when SNR is varied, the correlation between predictors is set to 0.8 times a random noise uniformly distributed on the interval [0.2, 0.4]. When the correlation ρ is varied, the SNR is fixed to 0.5.
It can be seen that as the SNR decreases or the correlation between the predictors increases, the estimation problem becomes harder, and the performance of all approaches gets worse. In general the Wasserstein ℓ 2 achieves the best performance in terms of all four metrics. Specifically, • It is better than the ℓ 2 -regularized LAD, which assumes an infinite transportation cost along y.
• It is better than the Wasserstein ℓ ∞ and ℓ 1 -regularized LAD which use the ℓ 1regularizer. • It is better than the approaches that use the SR loss function.
Empirically we have found out that in most cases, the approaches that use the AD loss, including the ℓ 1 -and ℓ 2 -regularized LAD, and the Wasserstein ℓ ∞ formulation, drive all the coordinates of β to zero, due to the relatively small magnitude of the AD loss compared to the norm of the coefficient, so that the regularizer dominates the solution. The approaches that use the SR loss, e.g., ridge regression and EN, do not exhibit such a problem, since the squared residuals weaken the dominance of the regularization term.
Overall the ℓ 2 -regularizer outperforms the ℓ 1 -regularizer, since the true regression coefficient is dense, which implies that a proper distance metric on the (x, y) space should take into account all the coordinates. From the perspective of the Wasserstein DRO framework, the ℓ 1 -regularizer corresponds to an ∥ · ∥ ∞ -based distance metric on the (x, y) space that only picks out the most influential coordinate to determine the closeness between data points, which in our case is not reasonable since every coordinate plays a role (reflected in the dense β*). In contrast, if β* is sparse, using the ∥ · ∥ ∞ as a distance metric on (x, y) is more appropriate. A more detailed discussion of this will be presented in Sections 4.3 and 4.4.

Dense β * , outliers only in x
In this subsection we will experiment with the same β* as in Section 4.1, but with perturbations only in x, i.e., for a given x of the outlier, the corresponding y value is drawn in the same way as the clean samples. Our goal is to investigate the performance of the Wasserstein formulation when the response y is not subjected to large perturbations. The motivation for introducing the AD loss in the Wasserstein formulation is to hedge against large residuals, as illustrated in Fig. 1. We are interested in comparing the AD and SR loss functions when the residuals have moderate magnitudes.
Interestingly, we have observed that although the ℓ 1 -and ℓ 2 -regularized LAD, as well as the Wasserstein ℓ ∞ formulation, exhibit unsatisfactory performance, the Wasserstein ℓ 2 , which shares the same loss function with them, is able to achieve a comparable performance with the best among all -EN and ridge regression (see Figs. 5 and 6). Notably, the ℓ 2 -regularized LAD, which is just slightly different from our Wasserstein ℓ 2 formulation, shows a much worse performance. This is because the ℓ 2 -regularized LAD implicitly assumes an infinite transportation cost along y, which gives zero tolerance to the variation in the response. For example, given two data points (x 1 , y 1 ) and (x 2 , y 2 ), as long as y 1 ≠ y 2 , the distance between them is infinity. Therefore, a reasonable amount of fluctuation, caused by the intrinsic randomness of y, would be overly exaggerated by the underlying metric used by the ℓ 2 -regularized LAD. In contrast, our Wasserstein approach uses a proper notion of norm to evaluate the distance in the (x, y) space and is able to effectively distinguish abnormally high variations from moderate, acceptable noise.
It is also worth noting that the formulations with the AD loss, e.g., ℓ 2 -and ℓ 1 -regularized LAD, and the Wasserstein ℓ ∞ , perform worse than the approaches with the SR loss.
penalizing the extended coefficient vector (−β, 1) seems to make up, making the Wasserstein ℓ 2 a competitive method even when the perturbations appear only in x.

Sparse β*, outliers in both x and y
In this subsection we will experiment with a sparse β*. The intercept is set to β 0 * = 3, and the coefficients for the 20 predictors are set to β* = (0.05, 0,0.006, 0, −0.007, 0, 0.008, 0, … , 0). The adversarial perturbations are present in both x and y. Specifically, the distribution of outliers is characterized by: 1.
Our goal is to study the impact of the sparsity of β* on the choice of the norm space for the Wasserstein metric. We know that the ℓ 1 -regularizer works better than the ℓ 2 -regularizer for sparse data, which has been validated by our results in Figs. 7 and 8. We will see that the Wasserstein ℓ ∞ formulation significantly outperforms the Wasserstein ℓ 2 . An intuitively appealing interpretation for the sparsity inducing property of the ℓ 1 -regularizer is made available by the Wasserstein DRO framework, which we explain as follows. The sparse regression coefficient β* implies that only a few predictors are relevant to the regression model, and thus when measuring the distance in the (x, y) space, we need a metric that only extracts the subset of relevant predictors. The ‖ · ‖ ∞ , which takes only the most influential coordinate of its argument, roughly serves this purpose. Compared to the ‖ · ‖ 2 which takes into account all the coordinates, most of which are redundant due to the sparsity assumption, ‖ · ‖ ∞ results in a better performance, and hence, the Wasserstein ℓ ∞ formulation that stems from the ‖ · ‖ ∞ distance metric on (x, y) and induces the ℓ 1 -regularizer is expected to outperform others.
We note that the ℓ 1 -regularized LAD achieves similar performance to ours, since replacing ‖β‖ 1 by ‖(−β, 1)‖ 1 only adds a constant term to the objective function. The generalization performance (mean MSE) of the AD loss-based formulations is consistently better than those with the SR loss, since the AD loss is less affected by large perturbations in y. Also note that choosing a wrong norm for the Wasserstein metric, e.g., the Wasserstein ℓ 2 , could lead to an enormous estimation error, whereas with a right norm space, we are guaranteed to outperform all others. Even when the SNR is very low, our performance is at least as good as the null estimator (see Fig. 7). Although EN and LASSO achieve similar performance to ours for moderate SNR values, they have a chance of performing even worse than the null estimator when there is little signal/information to learn from.

Sparse β*, outliers only in x
In this subsection, we will use the same sparse coefficient as in Section 4.3, but the perturbations are present only in x. Specifically, for outliers, their predictors and responses are drawn from the following distributions: 1.
Not surprisingly, the Wasserstein ℓ ∞ and the ℓ 1 -regularized LAD achieve the best performance. Notice that in Section 4.3, where perturbations appear in both x and y, the AD loss-based formulations have smaller generalization and estimation errors than the SR lossbased formulations. When we reduce the variation in y, the SR loss seems superior to the AD loss, if we restrict attention to the improperly regularized (ℓ 2 -regularizer) formulations (see Fig. 9). For the ℓ 1 -regularized formulations, our Wasserstein ℓ ∞ formulation, as well as the ℓ 1 -regularized LAD, is comparable with the EN and LASSO. Moreover, when there is little information to utilize (low SNR), EN and LASSO are worse than the null estimator, whereas our performance is at least as good as the null estimator.
We summarize below our main findings from all sets of experiments we have presented:

1.
When a proper norm space is selected for the Wasserstein metric, the Wasserstein DRO formulation outperforms all others in terms of the generalization and estimation qualities.

2.
Penalizing the extended regression coefficient (−β, 1) implicitly assumes a more reasonable distance metric on (x, y) and thus leads to a better performance.

3.
The AD loss is remarkably superior to the SR loss when there is large variation in the response y.

4.
The Wasserstein DRO formulation shows a more stable estimation performance than others when the correlation between predictors is varied.

An outlier detection example
As an application, we consider an unlabeled two-class classification problem, where our goal is to identify the abnormal class of data points based on the predictor and response information using the Wasserstein formulation. We do not know a priori whether the samples are normal or abnormal, and thus classification models do not apply. The commonly used regression model for this type of problem is the M-estimation (Huber, 1964(Huber, , 1973, against which we will compare in terms of the outlier detection capability.
The data are generated in the same fashion as before. For clean samples, all predictors x 1 , … , x 30 come from a normal distribution with mean 7.5 and standard deviation 4.0.
The response is a linear function of the predictors with β 0 * = 0.3, β 1 * = ⋯ = β 30 * = 0.5, plus a Gaussian distributed noise term with zero mean and standard deviation σ. The outliers concentrate in a cloud that is randomly placed in the interior of the x-space. Specifically, their predictors are uniformly distributed on (u − 0.125, u + 0.125), where u is a uniform random variable on (7.5 − 3 × 4, 7.5 + 3 × 4). The response values of the outliers are at a δ R distance off the regression plane. y = β 0 * + β 1 * x 1 + ⋯ + β 30 * x 30 + δ R .
the Receiver Operating Characteristic (ROC) curve which plots the true positive rate against the false positive rate, and the related Area Under Curve (AUC).
Notice that all the regression methods under consideration only generate an estimated regression coefficient. The identification of outliers is based on the residual and estimated standard deviation of the noise. Specifically, where σ is the standard deviation of residuals in the entire training set. ROC curves are obtained through adjusting the threshold value.
The regularization parameters for Wasserstein DRO and regularized LAD are tuned using a separate validation set as done in previous sections. We would like to highlight a salient advantage of our approach reflected in its robustness w.r.t. the choice of ϵ. In Fig. 11 we plot the out-of-sample AUC as the radius ϵ (regularization parameter) varies, for the ℓ 2 -induced Wasserstein DRO and the ℓ 1 -regularized LAD. For the Wasserstein DRO curve, when ϵ is small, the Wasserstein ball contains the true distribution with low confidence and thus AUC is low. On the other hand, too large ϵ makes our solution overly conservative. Note that the robustness of our approach, indicated by the flatness of the Wasserstein DRO curve, constitutes another advantage, whereas the performance of LAD dramatically deteriorates once the regularizer deviates from the optimum. Moreover, the maximal achievable AUC for Wasserstein DRO is significantly higher than LAD.
In Fig. 12 we show the ROC curves for different approaches, where q represents the percentage of outliers, and δ R the outlying distance along y. We see that the Wasserstein DRO formulation consistently outperforms all other approaches, with its ROC curve lying well above others. In general, all approaches have better performance when the percentage of outliers is lower, and the outlying distance is larger. The approaches that use the AD loss function (e.g., Wasserstein DRO and regularized LAD) tend to outperform those that adopt the SR loss (e.g., M-estimation which uses a variant of the SR loss). The superiority of our formulation could be attributed to the AD loss function, and the distributional robustness since we hedge against a family of plausible distributions, including the true distribution with high confidence. By contrast, M-estimation adopts an Iteratively Reweighted Least Squares (IRLS) procedure which assigns weights to data points based on the residuals from previous iterations, and then solves a weighted least squares estimation problem. With such an approach, there is a chance of exaggerating the influence of outliers while downplaying the importance of clean observations, especially when the initial residuals are obtained through Ordinary Least Squares (OLS).

Conclusions
Wasserstein metric was utilized to construct the ambiguity set and a tractable reformulation was derived. It is worth noting that the linear law assumption does not necessarily limit the applicability of our model. In fact, by appropriately pre-processing the data, one can often find a roughly linear relationship between the response and transformed explanatory variables. Our Wasserstein formulation incorporates a class of models whose specific form depends on the norm space that the Wasserstein metric is defined on. We provide out-of-sample generalization guarantees, and bound the estimation bias of the general formulation.
Extensive numerical examples demonstrate the superiority of the Wasserstein formulation and shed light on the advantages of the ℓ 1 -loss, the implication of the regularizer, and the selection of the norm space for the Wasserstein metric. We also presented an outlier detection example as an application of this robust learning procedure. A remarkable advantage of our approach rests in its flexibility to adjust the form of the regularizer based on the characteristics of the data. We want to find the set of θ such that the optimal values of problems A and B are finite.
Then, Dual-A and Dual-B need to have non-empty feasible sets, which implies the following two conditions: ∃r B ≥ 0, s.t. βr B = θ + β.

A.3 Proof of Theorem 3.3 Proof
We use Theorem 8 in Bartlett and Mendelson (2002), setting the following correspondences with the notation used there: ℒ(x, y) = ϕ(x, y) = y − x′β . This yields the bound (14) on the expected loss. For Eq. (15), we apply Markov's inequality to obtain:

A.4 Proof of Corollary 3.4 Proof
The percentage difference requirement can be translated into: from which (16) can be easily derived. ■

A.5 Proof of Corollary 3.5 Proof
Based on Theorem 3.3, we just need the following inequality to hold: which is equivalent to: We cannot obtain a lower bound for N by directly solving (24)

A.6 Sub-Gaussian Random Variables and Gaussian Width Definition 1 (Sub-Gaussian random variable)
A random variable z is sub-Gaussian if the ψ 2 -norm defined below is finite, i.e., z ψ 2 ≜ sup q ≥ 1 E z q q < + ∞ .
An equivalent property for sub-Gaussian random variables is that their tail distribution decays as fast as a Gaussian, namely, for some constant C.
A random vector z ∈ ℝ m is sub-Gaussian if z′u is sub-Gaussian for any u ∈ ℝ m . The ψ 2norm of a vector z is defined as: where S m denotes the unit sphere in the m-dimensional Euclidean space. For the properties of sub-Gaussian random variables/vectors, please refer to the book by Vershynin (2017).

Definition 2 (Gaussian width)
For any set A ⊆ ℝ m , its Gaussian width is defined as: where g N(0, I) is an m-dimensional standard Gaussian random vector.
On the other hand, from the Cauchy-Schwarz inequality: Combining (27) and (28), we have: where the last step follows from the fact that β est − β true / β est − β true 2 ∈ A β* . ■
Then, for any f w ∈ ℱ, where we used Γ = E zz′ and the fact that w ∈ A Γ .
For any f w ∈ ℱ we have f w ψ 2 = z′Γ −1/2 w ψ 2 where the first inequality used Assumption F and the second inequality used Assumption G.

A.9 Proof of Lemma 3.8
We follow the proof of Lemma 4 in Chen and Banerjee (2016), adapted to our setting. We include all key steps for completeness.
Then, by the tail behavior of sub-Gaussian random variables (see Hoeffding bound, Thm. 2.6.2 in (Vershynin, 2017)), we have: for some positive constant C 01 .
The result follows. ■

A.12 Proof of the Result in Section 4
We will show that if the Wasserstein metric is defined by the following metric s c : s c (x, y) = (x, cy) ∞ , then as c → ∞, the corresponding Wasserstein DRO formulation becomes: which is the ℓ 1 -regularized LAD.
We will use the notation z ≜ (x, y). Based on the definition of dual norm, we are interested in solving the following optimization problem for β ∈ ℝ m : max z z′β s.t. Z w, ∞ ≤ 1.